function [AB, AL, AC, AR, AT, S, dx, dy, dv, dc, siga, nusig, sct] = ...
    twoDcoefMGrespL( dat, abdat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, BC, ...
    legord,gg)
global order
% this simplified coefficient generator produces an incident current
% on only the bottom face of an assumed four-way symmetric box.
% The incident current is a DLP of order legord and group gg

% determine legendre polynomial for incoming current
N = sum(xfm);    % max degree
k = 0:N;         % interval
maxord = order;  % maxorder to which we compute P's
P = zeros(N+1,N+1);
P(1,1:N+1) = 1; 
P(2,1:N+1) = 1 - (2*k)/N;
for j = 3:N+1 %otherwise N-1
    i = j-2;
    P(i+2,1:N+1) = ( (2*i+1)*(N-2*k).*P(i+1,:) - i*(N+i+1)*P(i,:) ) ./ ...
        ((i+1)*(N-i)) ;
end
SS = zeros(N+1,numg);
SS(:,gg) = P(legord+1,:)';  % this is the boundary source.

% define data components
dc = zeros(numm,numg); sa=dc; ab=dc;
sc = zeros(numm,numg,numg);
for m = 1:numm
    for g = 1:numg %g+(m-1)*g
       dc(m,g) = dat(numg*(m-1)+g,1); % diffusion coefficient
       ab(m,g) = abdat(numg*(m-1)+g,1); % absorption xsec
       sa(m,g) = dat(numg*(m-1)+g,2); % removal cross-section
       ns(m,g) = dat(numg*(m-1)+g,3); % fission cross-section
       for gg = 1:numg
           sc(m,g,:) = dat(g+(m-1)*g,4:end);    
       end
   end
end

% number of fine meshes in x and y coordinates
N = sum(xfm); 
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm);
% compute dx and dy vectors
dx = zeros(N,1);
dy = zeros(M,1);
dv = zeros( (N+1)*(M+1), 1 );
j = 0;
for i = 1:CN
    dx( (j+1):(j+xfm(i)) ) = (xcm(i+1)-xcm(i))/xfm(i);
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    dy( (j+1):(j+yfm(i)) ) = (ycm(i+1)-ycm(i))/yfm(i);
    j = sum(yfm(1:i));
end

lenk = (N+1)*(M+1);
AC = zeros( lenk, numg );
AL = zeros( lenk-1, numg );
AR = AL;
AB = zeros( lenk-N-1, numg );
AT = AB;
S = zeros( (N+1)*(M+1), numg ); nusig=S; siga=S;
sct = zeros( (N+1)*(M+1), numg, numg );

for g = 1:numg % group loop

% interior coefficients
for i = 2:N 
    for j = 2:M
        for ci = 1:CN
            if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
        end
        for ci = 1:CM
            if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
        end
        for ci = 1:CN
            if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
        end
        for ci = 1:CM
            if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
        end        
        mtt  = mt(mtx,mty);   mtu  = mt(mtx,mtyu);
        mto  = mt(mtxo,mty);  mtb  = mt(mtxo,mtyu);
        
        k = i+(j-1)*( N+1 ); % maps i,j to k per duderstadt
       
        alft = -0.5*( dc(mtt,g)*dy(j-1) + dc(mtu,g)*dy(j) ) / dx(i-1);
        argt = -0.5*( dc(mto,g)*dy(j-1) + dc(mtb,g)*dy(j) ) / dx(i);
        abot = -0.5*( dc(mtt,g)*dx(i-1) + dc(mto,g)*dx(i) ) / dy(j-1);
        atop = -0.5*( dc(mtu,g)*dx(i-1) + dc(mtb,g)*dx(i) ) / dy(j);    
        a    = 0.25*( dx(i-1)*dy(j-1)*sa(mtt,g) + dx(i)*dy(j-1)*sa(mto,g) + ...
               dx(i-1)*dy(j)*sa(mtu,g) + dx(i)*dy(j)*sa(mtb,g) ) - (alft+argt+abot+atop);
        AL( k-1, g )   = alft;  
        AR( k, g )     = argt;   
        AB( k-N-1, g ) = abot;    
        AT( k, g )     = atop;
        AC( k, g )     = a;
        sct(k,g,:) =  0.25*( dx(i-1)*dy(j-1)*sc(mtt,g,:) + dx(i)*dy(j-1)*sc(mto,g,:) + ...
            dx(i-1)*dy(j)*sc(mtu,g,:) + dx(i)*dy(j)*sc(mtb,g,:) );
        nusig( k, g ) = 0.25*( dx(i-1)*dy(j-1)*ns(mtt,g) + dx(i)*dy(j-1)*ns(mto,g) + ...
            dx(i-1)*dy(j)*ns(mtu,g) + dx(i)*dy(j)*ns(mtb,g) );  
        siga( k, g )  = 0.25*( dx(i-1)*dy(j-1)*ab(mtt,g) + dx(i)*dy(j-1)*ab(mto,g) + ...
            dx(i-1)*dy(j)*ab(mtu,g) + dx(i)*dy(j)*ab(mtb,g) ); 
        dv( k, 1 ) = 0.25*( dx(i-1)*dy(j-1) + dx(i)*dy(j-1) + ...
            dx(i-1)*dy(j) + dx(i)*dy(j) );        
    end
end

% left and right edges
for j = 2:M
    for ci = 1:CN
        if 2 <= sum(xfm(1:ci)), mtx = ci; break; end
    end
    for ci = 1:CM
        if j-1 <= sum(yfm(1:ci)), mty = ci; break; end
    end
    for ci = 1:CN
        if N-1  <= sum(xfm(1:ci)), mtxo = ci; break; end
    end
    for ci = 1:CM
        if j  <= sum(yfm(1:ci)), mtyu = ci; break; end
    end       
    % left side
    i = 1;
    k = i+(j-1)*( N+1 );
    dc1 = 0.5*( dc(mt(1,mty),g) + dc(mt(1,mtyu),g) );
    dc2 = ( dx(i)*dy(j-1)*dc(mt(1,mty),g) + dx(i+1)*dy(j-1)*dc(mt(mtx,mty),g) +...
            dx(i)*dy(j)*dc(mt(1,mtyu),g)  + dx(i+1)*dy(j)*dc(mt(1,mtyu),g) ) / ...
          ( dx(i)*dy(j-1) + dx(i+1)*dy(j-1) + dx(i)*dy(j) + dx(i+1)*dy(j) );
    AR( k, g ) = -0.5 * dc2 / dx(i);    
    AC( k, g ) = 0.25 + 0.5 / dx(i) * dc1;
    % right side
    i = N+1;
    k = i+(j-1)*( N+1 );
    dc1 = 0.5*( dc(mt(end,mty),g) + dc(mt(end,mtyu),g) );
    dc2 = ( dx(N)*dy(j-1)*dc(mt(end,mty),g) + dx(N-1)*dy(j-1)*dc(mt(mtxo,mty),g)  +...
            dx(N)*dy(j)*dc(mt(end,mtyu),g)  + dx(N-1)*dy(j)*dc(mt(mtxo,mtyu),g) ) / ...
          ( dx(N-1)*dy(j-1) + dx(N)*dy(j-1) + dx(N-1)*dy(j) + dx(N)*dy(j) );
    AL( k-1, g ) = -0.5 * dc2 / dx(N);
    AC( k, g ) = 0.25 + 0.5 / dx(N) * dc1;                                     
end

% top and bottom edges
for i = 2:N
    for ci = 1:CN
        if i-1 <= sum(xfm(1:ci)), mtx = ci; break; end
    end
    for ci = 1:CM
        if 2 <= sum(yfm(1:ci)), mty = ci; break; end
    end
    for ci = 1:CN
        if i  <= sum(xfm(1:ci)), mtxo = ci; break; end
    end
    for ci = 1:CM
        if M-1  <= sum(yfm(1:ci)), mtyu = ci; break; end
    end 
    % bottom
    j = 1;
    k = i+(j-1)*( N+1 );  
    AT( k, g )     = -0.25 * ( dc(mt(mtx,1),g) + dc(mt(mtxo,1),g) ) / dy(1);   
    AC( k, g )     = 0.00 + 0.5 / dy(1) * ...
        ( dx(i-1)*dy(j)  *dc(mt(mtx,1),g)    + dx(i)*dy(j)  *dc(mt(mtxo,1),g) + ...
        dx(i-1)*dy(j+1)*dc(mt(mtx,mty),g)  + dx(i)*dy(j+1)*dc(mt(mtxo,mty),g) ) / ...
        ( dx(i-1)*dy(j) + dx(i)*dy(j) + dx(i-1)*dy(j+1) + dx(i)*dy(j+1) );
    S( k, g ) = SS(i,g); 
    % top
    j = M+1;
    k = i+(j-1)*( N+1 );
    AB( k-N-1, g ) = -0.25 * ( dc(mt(mtx,end),g) + dc(mt(mtxo,end),g) ) / dy(M);    
    AC( k, g )     = 0.25 + 0.5 / dy(M) * ...
        ( dx(i-1)*dy(M-1)*dc(mt(mtx,mtyu),g) + dx(i)*dy(M-1)*dc(mt(mtxo,mtyu),g) + ...
          dx(i-1)*dy(M)  *dc(mt(mtx,end),g)  + dx(i)*dy(M)  *dc(mt(mtxo,end),g) ) / ...
        ( dx(i-1)*dy(M-1) + dx(i-1)*dy(M) + dx(i)*dy(M-1) + dx(i)*dy(M) );                                   
end

% finally, the corners
for ci = 1:CN
    if i-1 <= sum(xfm(1:ci)), mtxL = ci; break; end  % picks out material for xfm(2)
end
for ci = 1:CM
    if 2 <= sum(yfm(1:ci)), mtyB = ci; break; end % picks out material for yfm(2)
end
for ci = 1:CN
    if i  <= sum(xfm(1:ci)), mtxR = ci; break; end % picks out material for xfm(end-1)
end
for ci = 1:CM
    if M-1  <= sum(yfm(1:ci)), mtyT = ci; break; end % picks out material for yfm(end-1)
end 
% BOTTOM LEFT
i=1; j=1; k = i+(j-1)*( N+1 );      
AR( k, g )     = -0.25 / dx(1) * ( dc(mt(1,1),g) + dc(mt(mtxL,1),g) ); 
AT( k, g )     = -0.25 / dy(1) * ( dc(mt(1,1),g) + dc(mt(1,mtyB),g) );
AC( k, g )     = 0.0 + 0.5*( dc(mt(1,1),g)/dx(1)+dc(mt(1,1),g)/dy(1) );
S( k, g )      = SS(i,g);   
% TOP LEFT
i=1; j=M+1; k = i+(j-1)*( N+1 );
AR( k, g )     = -0.25 / dx(1) * ( dc(mt(1,end),g) + dc(mt(mtxL,end),g) );  
AB( k-N-1, g ) = -0.25 / dy(end) * ( dc(mt(1,end),g) + dc(mt(1,mtyT),g) ); 
AC( k, g )     = 0.25 + 0.5*( dc(mt(1,end),g)/dx(1) + dc(mt(1,end),g)/dy(end) );       
% BOTTOM RIGHT
i=N+1; j=1;  k = i+(j-1)*( N+1 );     
AL( k-1, g )   = -0.25 / dx(end) * ( dc(mt(end,1),g) + dc(mt(mtxR,1),g) ) ;  
AT( k, g )     = -0.25 / dy(1) * ( dc(mt(end,1),g) + dc(mt(end,mtyB),g) ); 
AC( k, g )     = 0.0 + 0.5*( dc(mt(end,1),g)/dx(end)+dc(mt(end,1),g)/dy(1) );   
S( k, g )      = SS(i,g);   
% TOP RIGHT
i=N+1; j=M+1;  k = i+(j-1)*( N+1 );   
AL( k-1, g )   = -0.25 / dx(end) * ( dc(mt(end,end),g) + dc(mt(mtxR,end),g) );  
AB( k-N-1, g ) = -0.25 / dy(end) * ( dc(mt(end,end),g) + dc(mt(end,mtyT),g) );
AC( k, g )     = 0.25 + 0.5*( dc(mt(end,end),g)/dx(end) + dc(mt(end,end),g)/dy(end) );          

end % group loop
